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Rigorous computer simulations of propagating electromagnetic fields have become an important tool for optical 
metrology and design of nanostructured optical components. A vectorial finite element method (FEM) is a good 
choice for an accurate modeling of complicated geometrical features. However, from a numerical point of view 
solving the arising system of linear equations is very demanding even for medium sized 3D domains. In numerics, 
a domain decomposition method is a commonly used strategy to overcome this problem. Within this approach 
the overall computational domain is split up into smaller domains and interface conditions are used to assure 
continuity of the electromagnetic field. Unfortunately, standard implementations of the domain decomposition 
method as developed for electrostatic problems are not appropriate for wave propagation problems. In an earlier 
paper we therefore proposed a domain decomposition method adapted to electromagnetic field wave propagation 
problems. In this paper we apply this method to 3D mask simulation. 

Keywords: 3D EMF simulations, microlithography, adaptive high-order finite-element method, FEM, multiple 
scattering 



This paper adresses on the rigorous simulation of the electromagnetic field within a cutout of a microlithography 
mask or semiconductor wafer. It might be confusing that the notion "domain decomposition method" appears 
in the numerical literature about Maxwell's equations and in the mask simulation community with - at a first 
glance - different meanings. In the mask simulation community, Adam and Neureuther introduced the domain 
decomposition method in 200 1. 1 As an example, in the approach by Adam and Neureuther a phase shift line 
mask is split up into two isolated lines, c.f. Figure[T] The transmission through each line is computed separately. 
Afterwards the computed fields are merged in an appropriate way. Diffraction effects of the double slit are 
neglected in a first step but can also be approximated by adding further cross talk terms. In the numerics 
community the domain decomposition method denotes a strategy for solving large scale partial differential 
equation problems and goes back to Schwarz more than 100 years ago. 2 This approach shares the key idea 
by Adam and Neureuther to split the computational problem into smaller, more feasible subproblems. But in 
contrast to the method by Adam and Neureuther the synthesis step (merging the fields on the sub-domains) 
relies on an iterative numerical method and not on physical intuition. This way the method remains rigorous: 
Starting from a rigorous discretization of Maxwell's equations in the entire domain one gets a large scale linear 
system of equations. The aim is to iteratively solve this large scale linear system. As long as the iterative process 
converges this yields the rigorous solution. Although the domain decomposition method is well established and 
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ABSTRACT 



1. INTRODUCTION 





Figure 1. Domain decomposition method by Adam and Neureuther (not applied in this paper). One uses a superposition 
of the fields Ei, E2, and E3 to approximate the field E. Although this approach yields satisfactory results for many 
applications in mask design it introduces errors inversely proportional to the distance of the features. 




Figure 2. Domain decomposition method as discussed in this paper (horizontal setting). For each sub-domain transparent 
boundary conditions are imposed. Starting with the first sub-domain the field E^ ' and the scattered field E^ c to the 
incoming wave Ei nc are computed. Afterwards the second domain is updated, where the incoming field is now the sum 
of the originally incident field Ei nc and the scattered field E^ c . This process is iteratively repeated. Observe that the 
zeroth iterate is exactly the method by Adam and Neureuther. 

well analyzed for electrostatic and eddy current problems the application of domain decomposition techniques 
to high frequency Maxwell's equations is still an active research area in the numerics community. Especially, 
for wave propagation problems it is necessary to impose proper radiation boundary conditions at the coupling 
boundary of the adjacent sub-domains. 3,4 We will detail this issue in the next section. Figure [2] and Figured] 
give a rough sketch of the algorithm. In both examples the computational domain is split into two isolated 
structures. Starting with the first sub-domain the field and the scattered field E^ to the incoming wave 
Ej llc are computed. Afterwards the second domain is updated, where the incoming field is now the sum of the 
originally incident field Ej nc and the scattered field Ej°g C . This process is iteratively repeated. This approach 
resembles the multiple scattering method. 5 Also observe that the zeroth iterate is exactly the method by Adam 
and Neureuther. 

The paper is structured as follows. Section[2]is devoted to the modeling of Maxwell's equations and of the domain 
decomposition iteration. We further summarize various numerical issues such as the definition of transparent 
boundary conditions and casting Maxwell's equations into weak form needed for the finite element method. The 
derivation of a weak form for Maxwell's equations on structured unbounded domains was already discussed in 
our papers, 6 ' 7 but also see the recent paper by Wei et al. 8 In section [3] we apply the domain decomposition 
method to a medium sized cutout of a mask. The paper ends with some concluding remarks. 

2. MODELING 

2.1 Maxwell's equations 

Starting from Maxwell's equations in a medium without sources and free currents and assuming time-harmonic 
dependence with angular frequency lu > the electric and magnetic fields 

E(x, y, z, t) = E(x, y, z)e~^\ H(x, y, z, t) = U(x, y, z)e~ l ^\ 

must satisfy 

V x E = iu/jH, V • eE = 0, 

V x H = -iweE, V • ^tH = 0. 
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Figure 3. Domain decomposition method as discussed in this paper (vertical setting). The notation is given in Figure 

Here e denotes the permittivity tensor and fj, denotes the permeability tensor of the materials. In the following we 
drop the wiggles, so that E — > E, H — ► H. From the equations above we then may derive (by direct substitution) 
the second order equation for the electric field 

V x /i _1 V x E - lu 2 £~E = 0, 
V • eE = 0. 

A similar equation holds true for the magnetic field - one only need to replace E by H and interchange e and 
\l. Observe that any solution to the first equation also meets the divergence condition (second equation). We 
therefore drop the second equation in the following. 

For the sake of a simpler and more compact notation we rewrite these equations in differential form, 

d 1 fi~ 1 d 1 e - uj 2 ee = 0. (1) 

A reader not familiar with this calculus may replace the exterior derivatives do, d\, g? 2 with classical differential 
operators, do — > V, di — > Vx and d 2 — > V-. Here, the electric field appears as a differential 1-form, e = 
e x dx + e y dy + e z dz, whereas the material tensors act - from a more mathematical point of view - as operators 

e, fi : Alt 1 -> Alt 2 . 

In the following we drop the sub- indices for the exterior derivatives do, d\, d 2 . 
2.2 Scattering off an isolated structure 

We now deal with light scattering off an isolated structure as depicted in Figured] The block fl — [—a, a] x 
[— b, b] x [— c, c] is the computational domain. For simplicity we assume that the computational domain is aligned 
along the coordinate axes and is embedded into a layered media. The role of the hyperplane T will be explained 
later. The incident field Ej nc is a solution to Maxwell's equations in the layered media. Typically Ej nc is a plane 
wave traveling in the z-direction scattered at the multi-layer stack. Outside the computational domain the total 
field is a sum of the incident field and the scattered field E sc , which also is a solution to Maxwell's equations in 
the layered media outside the computational domain. The interior and exterior fields are linked by the following 
coupled system, 

d/j~ 1 dej nt — u! 2 eej nt = 0, on ft (2a) 
d[i~ 1 de sc - uj 2 ee sc = 0, on R 3 \ ft (2b) 
eint - e sc = e inc , on <9f2 (2c) 
A*int lrfe mt - MtTxt^sc = AWtCtew, on dtt. (2d) 



The first two equations are Maxwell's equations for the interior and the scattered field. Equation ([2"c|) and 
equation (|2d| assure the tangential continuity of the electric field and the magnetic field respectively. 



Figure 4. Sketch of a 3D computational domain. The structure is embedded into a multi-layer stack (mask or wafer 
blank). The incoming field Ei nc is a solution to Maxwell's equation in the layered media. Typically a plane wave traveling 
in the z-direction and scattered at the multi-layer stack is used. Periodic problems can also be dealt with but are not 
specially detailed in this paper. The plane V is used in the domain decomposition process. Within this context the 
sketched domain is only a sub-domain of the entire domain. F separates two regions which are merged in the domain 
decomposition iteration. 



2.3 Domain decomposition iteration 

In the domain decomposition method the entire computational domain is split into n smaller boxes fli. Initially, 
one starts with one sub-domain and solves the coupled interior-exterior system ([2]). When updating the subse- 
quent domains we account for the already computed scattered field of the adjacent domains. 
Let us first regard the instructive example as in Figure [5] There the entire domain is split into three sub-domains 
aligned in a linear chain. Assume that initially only the scattered field on domain O3 is computed. Next the 
middle domain Q 2 is updated. We solve system © where the incident field is now the sum of the original light 
source field and the scattered field E SCj 3 from domain Since the domain SI2 does not contain any scatterer 
the scattered field E sc ,2 from this domain is zero. To update domain f2i now we have to gather all waves which 
travel through interface r^. Obviously, this field consists of the field scattered from domain travelled through 
^2 and the original source light from the illumination system. The key issue here is, that it is not sufficient to 
compute the scattered field for (which is zero here). We also need to compute the field travelling from one 
parallel interface to the other. 

We now deal with light scattering off an isolated sub-domain as depicted in Figure |U We introduce the infinite 
hyperplanes Tij which separates the two sub-domains Qi and flj with normal mj directed from Qi to fl,-, see 
hyperplane T in Figure Q] We define the restriction of the scattered field e SCji onto the hyperplane 

These data are stored in the domain decomposition iteration, so that in the fcth iteration data r . . are given. 
For the update step 

_(*) _^ e (*+i) 




Figure 5. Domain decomposition of the entire domain in three sub-domains aligned in a line. As an instructive example, 
the middle domain does not contain a scatterer. The field E SCj 3 scattered from domain Q3 enters the middle domain 
across the infinite interface F2,3. In the domain decomposition process the field travels through the middle domain and 
hits domain f2i across the interface Ti 2- 



we solve the following system on fii which accounts for propagation through the sub-domains as well as scattering 
within the sub-domains 

d/i~ 1 dei Dt — uj 2 eei nt — 0, on Q (3a) 
dn~ 1 de sc - uj 2 ee sc = 0, on R 3 \(QU T it .) (3b) 

eint - e sc = e inc + e^.. , on dfl n T ltj (3c) 

Mhrtde int - Vcxt de sc = Mext^inc + M^xt^sc^.i > on dfl f] T itj (3d) 

e sc _-e sc , + = e^ ., on T it j \ dfl (3e) 
VcAde sc - - fJ-^ t de SCt+ = /x^de^. ., on T ij3 \ dfl. (3f) 

With equation (|3b[) the so defined scattered field e sc meets Maxwell's equations in the exterior domain but may 
jump across the hyperplanes r^. For each hyperplane Tij the jump is defined in equations (f3"c|) . (fSejl and is 
equal to the scattered field of the corresponding adjacent domain. On r,j the quantity e SCi+ is the field limit in 
the riij-direction and e sc ,_ is the field limit in the opposite direction. As the update step we now define 

e [k+1) - (e ^ 

From equations (|3"c|) and (|3d|) one shows that under convergence of the iteration the field continuity across 
Yi.j fl dfl is satisfied. Hence the so merged field satisfies Maxwell's equations on the entire domain. 

2.4 Weak formulation 

To apply the finite element method we need to transform the linear systems used in the previous section to 
a so called weak form. The incorporation of the various jump conditions in system ^ is very technical, but 
fits seamless into the finite element framework. We exemplify this only for the simpler system of equations (0). 
The treatment of the system ([3]) will be commented at the end of sub-section 12.51 fn order to derive a weak 
formulation we define the following function space on the domain fl = R 3 

#i oc (curl) = {e G Alt 1 | (e x ,e y , e z ) e (^L) 3 , V x (e x ,e y , e 2 ) T 6 (L^ oc ) 3 } . 

The weak form of equations now reads 



/ (fi~ 1 de A dv- uj 2 (ee) A v) = 

R 3 



(4) 



for all v € 7?i oc (curl) with compact support. We now look for a variational formulation, where the incoming 
field only appears on the right hand side. To do that let us denote the computational domain by Q and split 
Maxwell's equations ([U into an interior and exterior part, 



/ (fi 1 de int A dv - Lu 2 (ee int ) A v) + / (/i 1 de sc A dv - uj 2 (ee sc ) A v) = 
Jn Jn 3 \n 

- / (// _1 <fei nc A dv — cu 2 (ee inc ) A v) 
Jr 3 \q 



(Sint e inc)|^0 — ( e sc)|9Q 



Applying a partial integration on the right hand yields 

/ (/i _1 de int Arfv-o; 2 (eei n t) A v) + / (// _1 cfe sc A dv — w 2 (ee sc ) A v) = / (^ 1 de inc A v) (5a) 
Jo jR 3 \n Jan 

(eint — e inc)|an = ( e sc)|an- (5b) 

On the left hand side we find the quantities of interest we want to compute, namely the total field in the interior 
ei n t and the scattered field e sc in the exterior. The second equation (|5bj) merges the scattered and the interior 
field. After introducing transparent boundary conditions we will explain how to incorporate this field data 
matching condition (|5bp into a variational formulation 

2.5 Transparent boundary condition (PML) 

So far, in all considerations the various scattering problems were posed on the entire domain R 3 and are therefore 
numerically not feasible. This is overcome by using transparent boundary conditions. We use the perfectly 
matched layer method introduced by Berenger. 9 This method exploits the analytic continuation properties of 
the scattered field in the exterior domain. In a nutshell using an appropriate complex continuation, the scattered 
field is transformed to an exponentially decaying field without affecting the matching condition with the field in 
the interior domain. In an earlier paper we proposed an extremely efficient adaptive PML method which also 
copes with non cubic domains and various kinds of inhomogeneous exterior domains. 6,7 For simplicity, in this 
paper we restrict ourselves to cubic domains embedded into an layered media. Let us regard the computational 
domain ft = [—a, a] x [—6, b] x [— c, c] in Figure 2J To derive the PML equation we use different coordinate 
stretchings in each coordinate direction, 10 e.g. 

a + j(x — a), x > a 

x, \x\ < a . 

—a + j(x + a), x < — a 

The definitions for the y and z-dircctions are accordingly. For a complex coordinate stretching 7 is a complex 
number with 5R7 = 1.0 and 37 > 0.0. But firstly we consider 7 as a real number. Then stretching the coordinates 
is a simple coordinate transformation. In the differential form calculus, when switching the coordinates, the 
differential forms (field data) and the operators e and \x are transformed accordingly. Subscribing 7 to the 
transformed quantities equation §5§ yields 

/ (^^cfent A dv - uj 2 (ee int ) A v) + / (/i~ 1 tfe SC)7 A dv^j — w 2 (ee SC!7 ) A v^) = / (^ _1 de inc A v) 
Jo, JR 3 \a Jan 

( e int — e inc) |gQ = ( e sc,7 ) \qq ■ 

This equation is identical to (|5|) but only the transformed quantities are used. Since the test function v is chosen 
arbitrarily we can avoid using the transformed field and replace v 7 by v without changing the variational form. 
Now we switch to a complex coordinate stretching. Since the equation above holds true for any real 7 and due 
to the holomorphy of the scattered field e sc it is a simple matter of complex function theory that the above 
equations also hold true for 7 chosen complex. We now want to incorporate the matching condition on the 
boundary dfl into the variational form. At the boundary Oil the interior field ei nt and the scattered field e SCi7 



differ by the field values of the incident field. We therefore add a field with tangential data equal to ej nc on the 
boundary dfl and which has local support in the exterior domain R 3 \ fl. Since this field interpolates the data 
of the incident field on the boundary we denote it by Ie mc . In the finite element context I is the boundary 
interpolation operator. With the definition e = ej nt + e sc-7 +2e- mc we get the variational formulation 

/ (^deAdv-w^eEjAv) = / lnc [v] (6) 
Jn 

with 

/inc[v] = - / (fi~ 1 dle inc A dv - uj 2 (ele inc ) A v) + / (/i _1 de inc A v) . 
Jn x Jon 

Here, ilj denotes the finite support of Ie mc . Hence on the right hand side only data of the incident field around 
the boundary <9f2 of the computational domain are involved. Truncating the exterior domain this variational 
form is suitable for a finite element discretization. In the following we again drop the wiggles. 
As we mentioned above deriving the weak form for the system used in the domain decomposition iteration @ is 
more tedious. In this case interior boundary conditions on the infinite hyperplanes Tij must be treated. After 
defining the transparent boundary conditions by the PML method these jump conditions are imposed within the 
PML sponge layer. Since the effect of the artificial interior boundaries vanishes with convergence of the iterative 
process the transparent boundary conditions on these interfaces only require less accuracy. Even approximate 
transparent boundary conditions can be used. 

3. NUMERICAL EXAMPLE: PERIODIC 3D MASK 

In this section we apply the domain decomposition method to a 3D mask cutout given in Figure El The structure 
consists of MoSi-lines of height h — 65.4nm on a glass substrate. Here, the sidewall angle of the lines is 90deg, 
but other sidewall angle parameter can be treated without extra costs. Further material values are given in 
the description to Figure [51 We simulated the mask with high order finite elements of order 5 and 6. Using 
finite clement degree 5 (1.6 million unknowns) it was still possible to compute the solution with the direct 
method PARDISO. 11 A comparison of the domain decomposition method and the direct method shows a good 
agreement. After five iterations the domain decomposition method with 16 vertical sub-domains converged to 
the exact solution up to an error of 10 -5 in the complex amplitudes of the diffraction modes. The required 
memory was reduced by a factor of 3 (from 71GB to 25GB). On our computer it was not possible to compute 
the solution to finite element degree 6 by a direct sparse LU method due to limited memory resources. However, 
the convergence rate of the domain decomposition method was not affected when increasing the finite element 
degree from 5 to 6. This way we were able to solve the linear system with more than 3.3 million of unknowns. 
However with 50GB needed RAM the memory requirements are still demanding. To reduce further the memory 
demand it seems promising to split the domains also in the horizontal direction, which will be done in future 
work. In Figure[7|the computed near fields for an perpendicular incidence of x— and y— polarized light is plotted 
in a distance of lbnm above the structure. 

4. CONCLUSIONS 

We have proposed a new rigorous domain decomposition method for Maxwell's equations. This method allows 
for the reduction of needed computer resources. Especially it is possible to reduce significantly the amount of 
needed memory. Is has been shown that the method converges to the exact (discrete) finite element solution 
which corresponds to a large scale linear system of equations. The convergence rate behaviour is currently under 
investigation. 
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